function output = GPOPSEndpoint(input)
    % turn into minimization problem
    welfare = input.phase.integral(1) + input.phase.integral(2) + ...
              input.phase.integral(3) + input.phase.integral(4);
    output.objective      = - welfare;

    [Y_B_M, Y_B_R, Y_B_C, Y_E, Y_S] = compute_capital_returns(input);
    
    % event constraint (equalize return to robots)
    output.eventgroup(1).event = [Y_B_M - Y_B_R, ...
                                  Y_B_M - Y_B_C];
    
    % labor market clearing
    inpscale = input.auxdata.settings.inpscale;    
    
    
    integral_M = input.phase.integral(6);
    integral_R = input.phase.integral(7);
    integral_C = input.phase.integral(8);
    
    inp.L_M =   inpscale.L_M .*   input.parameter(1);
    inp.L_R =   inpscale.L_R .*   input.parameter(2);
    inp.L_C =   inpscale.L_C .*   input.parameter(3);
    inp.K_B_M = inpscale.K_B_M .* input.parameter(4);
    inp.K_B_R = inpscale.K_B_R .* input.parameter(5);
    inp.K_B_C = inpscale.K_B_C .* input.parameter(6);
    inp.K_E =   inpscale.K_E .*   input.parameter(7);
    inp.K_S =   inpscale.K_S .*   input.parameter(8);
    
    prod_fun = input.auxdata.gpops_fun.prod_fun;
    
    Y_M = prod_fun.Y_M(inp);
    Y_R = prod_fun.Y_R(inp);
    Y_C = prod_fun.Y_C(inp);    
    
    output.eventgroup(2).event = [integral_M - Y_M * inp.L_M, ...
                                  integral_R - Y_R * inp.L_R, ...
                                  integral_C - Y_C * inp.L_C];    
    
    
    settings = input.auxdata.settings;
    q_B = input.auxdata.gpops_fun.q_B;
    q_E = input.auxdata.gpops_fun.q_E;
    q_S = input.auxdata.gpops_fun.q_S;
    
    tau_B = input.auxdata.settings.tau_B;
    tau_E = input.auxdata.settings.tau_E;
    tau_S = input.auxdata.settings.tau_S;
    

    if settings.default_capital_tax
        % does not include default robot tax, only defaults for
        % equipment and structures
        output.eventgroup(1).event = [output.eventgroup(1).event, Y_E - (1+tau_E)*q_E,Y_S-(1+tau_S)*q_S];
    end
    
    if settings.default_robot_tax % option used for version with
                                  % default capital taxes
        output.eventgroup(1).event = [output.eventgroup(1).event ,Y_B_M - (1+tau_B)*q_B];
    end     
    
    if settings.equipment_equal_robot_tax
        output.eventgroup(1).event = [output.eventgroup(1).event, ...
                                      Y_B_M/q_B - 1 - (Y_E/q_E - 1)];
    end
    
end


function [Y_B_M, Y_B_R, Y_B_C, Y_E, Y_S] = compute_capital_returns(input)
    inpscale = input.auxdata.settings.inpscale;    
    
    inp.L_M =   inpscale.L_M .*   input.parameter(1);
    inp.L_R =   inpscale.L_R .*   input.parameter(2);
    inp.L_C =   inpscale.L_C .*   input.parameter(3);
    inp.K_B_M = inpscale.K_B_M .* input.parameter(4);
    inp.K_B_R = inpscale.K_B_R .* input.parameter(5);
    inp.K_B_C = inpscale.K_B_C .* input.parameter(6);
    inp.K_E =   inpscale.K_E .*   input.parameter(7);
    inp.K_S =   inpscale.K_S .*   input.parameter(8);
    
    prod_fun = input.auxdata.gpops_fun.prod_fun;
    
    Y_B_M = prod_fun.Y_B_M(inp);
    Y_B_R = prod_fun.Y_B_R(inp);
    Y_B_C = prod_fun.Y_B_C(inp);
    Y_E = prod_fun.Y_E(inp);
    Y_S = prod_fun.Y_S(inp);
end